library(tidyverse)
library(ggplot2)
theme_set(theme_bw())
library(xtable)
library(matrixStats)
library(descr)
library(readstata13)
library(readr)
library(jtools)
library(interactions)
library(modelsummary)
library(broom)
library(estimatr)
library(patchwork)
library(interflex)

set.seed(1234557)
data1 <- read_csv("PV_replication.csv")


data1<- data1 %>% 
  mutate(treatvar = as.factor(treatvar),
         treatdum = as.factor(treatdum))

data2<- data1 %>% 
  mutate(treatvar = as.numeric(treatvar),
         treatdum = as.numeric(treatdum))

data3<- data1 %>% 
  mutate(treatvar = as.factor(treatvar),
         treatdum = as.factor(treatdum),
         outcomecat = as.factor(outcome_1week))


##Figure 1##

coefficients_1week<-lm_robust(outcome_1week ~ treatdum, data=data1)$coefficients
coefficients_3days<-lm_robust(outcome ~ treatdum, data=data1)$coefficients

ses_1week<-(lm_robust(outcome_1week ~ treatdum, data=data1)$std.error)*100
ses_3days<-(lm_robust(outcome ~ treatdum, data=data1)$std.error)*100


tab0<-NA

tab0 <- data.frame(c(rep("t+3",1), rep("t+3",1), rep("t+7",1), rep("t+7",1)))
tab0$Time <- NA
tab0$Condition <- NA
tab0$mean <- NA
tab0$Cilower <- NA
tab0$Ciupper <- NA
colnames(tab0) <- c("Time", "Condition", "Mean", "Ci_lower", "Ci_upper")

tab0[1,3]<-coefficients_3days[1]*100
tab0[2,3]<-coefficients_3days[2]*100
tab0[3,3]<-coefficients_1week[1]*100
tab0[4,3]<-coefficients_1week[2]*100

tab0$Ci_upper<- c(coefficients_3days[1]*100 + 1.96*ses_3days[1], coefficients_3days[2]*100 + 1.96*ses_3days[2], coefficients_1week[1]*100 + 1.96*ses_1week[1], 
                  coefficients_1week[2]*100 + 1.96*ses_1week[2])

tab0$Ci_lower<- c(coefficients_3days[1]*100 - 1.96*ses_3days[1], coefficients_3days[2]*100 - 1.96*ses_3days[2], coefficients_1week[1]*100 - 1.96*ses_1week[1], 
                  coefficients_1week[2]*100 - 1.96*ses_1week[2])

tab0[1,2]<-"Control"
tab0[2,2]<-"Email"
tab0[3,2]<-"Control"
tab0[4,2]<-"Email"


ggplot(tab0, aes(x=Time, y=Mean, colour=Condition)) + 
  geom_errorbar(aes(ymin=Ci_lower, ymax=Ci_upper), width=.07) +  scale_colour_manual(values=c("grey40", "mediumvioletred", "grey40", "mediumvioletred")) +
  geom_line() + geom_point(aes(shape=Condition)) + geom_hline(yintercept=0, colour="grey", linetype = "dashed", size = 1) +
  theme_bw() + ggtitle("% writing to their MP after 3 and 7 days") + coord_cartesian(ylim = c(0, 4)) +  ylab("% writing to their MP") +  xlab("Days after treatment")
ggsave("Figure1.png", dpi=600)


##Figure 2##
model1  <- lm(outcome_1week ~ treatvar, data=data1)
msummary(model1, star=TRUE, output='latex')
data1$predicted<-predict(model1, data1)

effect_plot(model = model1, pred = treatvar, 
            cat.geom="point", cat.interval.geom="linerange", robust=TRUE)+
  labs(title = "Proportion of respondents who contacted their MP (95% CIs)")+
  ylab("Pr(Contacted their MP)")+
  xlab("Treatment")+
  scale_x_discrete(labels=c("Control" = "Control", "Email_noposition_nourgency" = "T1 \nno position & no urgency",
                            "Email_noposition_urgency" = "T2 \nno position & urgency", 
                            "Email_position_nourgency" = "T3 \nposition & no urgency ", 
                            "Email_position_urgency" = "T4 \nposition & urgency"))+
  geom_jitter(data = filter(data1, outcome_1week>0),
              aes(x=treatvar, y=predicted),
              size=1, alpha=.25, width=.4, height=.01, 
              pch=21, color="white", fill="mediumvioletred")+
  annotate(
    geom="text", x = 2.5, y = 0.015, size = 3, color = "grey57", 
    label = "Respondent who \ncontacted their MP")+
  annotate(
    geom = "curve", x = 2.5, y = .013, xend = 1.3, yend = .006, 
    size=.5, curvature = -.2, arrow = arrow(length = unit(2, "mm")), colour="grey57")

ggsave("Figure2.png", dpi=600)


##Figure 3##
datasub <- read_csv("WTT_data.csv")
subplot<- ggplot(datasub, aes(x = running, y = count, color=party)) +
  scale_color_manual(values = c("purple1", "blue2", "red2", "orange2", "green4"), 
                     labels = c("All MPs", "Conservative", "Labour", "Lib Dem", "Others"))+
  geom_smooth(method = "gam", show.legend = FALSE) +
  geom_jitter(data = datasub, size=1, alpha=.8) +
  geom_vline(xintercept = 0)+
  theme_bw()+
  ylim(0, 420)+
  annotate(
    geom="text", x = 20, y = 410, size = 3, color = "gray57", fontface=2,
    label = "People's Vote \ncampaign email")+
  annotate(
    geom = "curve", x = 20, y = 390, xend = 0, yend = 310, 
    curvature = -.2, arrow = arrow(length = unit(2, "mm")), colour="gray57")+
  labs(x = "Day relative to treatment", 
       y = "Daily MP contact \nvia WriteToThem.com",
       title="(No) evidence of substitution effects",
       color = "")+
  theme(legend.position = "bottom")+
  guides(color=guide_legend(nrow=1, byrow=TRUE))+
  theme(
    plot.title = element_text(hjust = 0, color = "gray57", size = 12, face="bold"),
    legend.text = element_text(hjust = 0, color = "gray57", size = 10, face="bold"),
    axis.title = element_text(hjust = .5, color = "gray57", size = 10, face="bold"))

ggsave("Figure3.png", dpi=600)

##Figure 4##

data1$donor_num=as.numeric(data1$donor)
data1$isdonor<-NA
data1$isdonor[data1$donor_num==1]<-1

data1$nodonor<-NA
data1$nodonor[data1$donor_num==0]<-1


sub_donor<-subset(data1, !is.na(isdonor))
sub_nodonor<-subset(data1, !is.na(nodonor))


full_coef<-lm_robust(outcome ~ position + urgency, data=data1)$coefficients
full_cilow<-lm_robust(outcome ~ position + urgency, data=data1)$conf.low
full_cihigh<-lm_robust(outcome ~ position + urgency, data=data1)$conf.high


donor_coef<-lm_robust(outcome ~ position + urgency, data=sub_donor)$coefficients
donor_cilow<-lm_robust(outcome ~ position + urgency, data=sub_donor)$conf.low
donor_cihigh<-lm_robust(outcome ~ position + urgency, data=sub_donor)$conf.high


nodonor_coef<-lm_robust(outcome ~ position + urgency, data=sub_nodonor)$coefficients
nodonor_cilow<-lm_robust(outcome ~ position + urgency, data=sub_nodonor)$conf.low
nodonor_cihigh<-lm_robust(outcome ~ position + urgency, data=sub_nodonor)$conf.high


tab3<-NA
tab3 <- data.frame(c(rep("1. Full sample",1), rep("2. Donor",1), rep("3. No donor",1)))
tab3$donor <- NA
tab3$ITT<- NA
tab3$Cilower <- NA
tab3$Ciupper <- NA
colnames(tab3) <- c("Donor", "ITT","Ci_lower", "Ci_upper")
tab3$ITT<- c(full_coef[2], donor_coef[2], nodonor_coef[2])*100
tab3$Ci_upper<- c(full_cihigh[2], donor_cihigh[2], nodonor_cihigh[2])*100
tab3$Ci_lower<- c(full_cilow[2], donor_cilow[2], nodonor_cilow[2])*100

graphdb <- ggplot(tab3,aes(x = Donor, y = ITT,ymin = Ci_lower, ymax = Ci_upper))
graphdb + scale_colour_manual(values=c("mediumvioletred", "mediumvioletred")) + geom_point(position=position_dodge(width=0.8) ,size = 3, color="mediumvioletred") + 
  geom_linerange(position=position_dodge(width=0.5), color="mediumvioletred", size =0.5) +
  geom_hline(yintercept=0, colour="grey", linetype = "dashed", size = 1) +
  theme_bw() + ggtitle("ITT of revealing position on writing to MP") + coord_cartesian(ylim = c(-1.5, 1.5)) +  ylab("ITT in %-points")
ggsave("Figure4.png", dpi=600)


##Figure 5##
#figure_touse variable from Hanretty (2017)#

model2  <- lm(outcome_1week ~ treatdum*figure_touse, data=data1)

interact1<- interact_plot(model = model2, pred =figure_touse, modx =treatdum, rug=TRUE,
                          interval=TRUE, colors="Rainbow", robust=TRUE, cluster="constituency",
                          modx.labels=c("0" = "Control", "1" = "Treatment"))+
  ylim(0, 0.05)+
  labs(title="Email receipt",
       y="Pr(Contact MP)", x="Constituency % leave vote")+
  theme(legend.position = c(0.4, 0.45),
        legend.title = element_blank())
ggsave("interaction1.png")


model3  <- lm(outcome_1week ~ treatvar*figure_touse, data=data1)
interact2<- interact_plot(model = model3, pred =figure_touse, modx =treatvar, rug=TRUE,
                          interval=TRUE, colors="Blues", robust=TRUE, cluster="constituency",
                          modx.labels=c("1" = "Control", "2" = "T1",
                                        "3" = "T2", "4" = "T3",
                                        "5" = "T4"))+
  ylim(0, 0.05)+
  labs(title="Heterogenous treatment messages",
       y="Pr(Contact MP)", x="Constituency % leave vote")+
  theme(legend.position = c(0.8, 0.35),
        legend.title = element_blank())
ggsave("interaction2.png")

interact3<- interact1 +interact2+
  plot_annotation(title = "Treatment effect moderated by constituency support for Brexit",
                  caption= "Model reporting constituency-clustered robust standard errors",
                  theme=theme(plot.title=element_text(face="bold", hjust=0)))

ggsave("Figure5.png", dpi=600)
